In [1]:
import plotly.graph_objects as go
import numpy as np

The velocity of a flow moving through a bend can be represented in polar coordinates: \begin{equation} \require{color}\vec{{\color[rgb]{0.059472,0.501943,0.998465}v}} = {\color[rgb]{0.059472,0.501943,0.998465}v}_{r}\hat{e}_{r} + {\color[rgb]{0.059472,0.501943,0.998465}v}_{\theta}\hat{e}_{\theta} \end{equation}

Because the flow is steady the streamlines will follow a circular path, therefore the radial component of velocity is zero, and the velocity will point in the angular direction

In [2]:
# initialize theta values for graph 
t1=np.linspace(-np.pi/8,np.pi/8,50)
t2=np.linspace(-np.pi/2, np.pi/2,50)
In [3]:
# plot duct and force diagram for small area of flow 
fig = go.Figure()
fig.add_scatter(x=2*np.cos(t2), y=2*np.sin(t2), fill='tozeroy', fillcolor='rgba(133, 163, 201, 0.5)', line_color='black')
fig.add_scatter(x=np.cos(t2), y=np.sin(t2), fill='tozeroy', fillcolor='white', line_color='black')

fig.add_scatter(x=1.25*np.cos(t1), y=1.25*np.sin(t1), line_color='black')
fig.add_scatter(x=1.75*np.cos(t1), y=1.75*np.sin(t1), fill='tonexty', line_color='black', fillcolor = 'steelblue')
fig.add_scatter(x=np.array([1.25,1.75])*np.cos(np.pi/8), y=np.array([1.25,1.75])*np.sin(np.pi/8), line_color='black', 
                mode='lines')
fig.add_scatter(x=np.array([1.25,1.75])*np.cos(-np.pi/8), y=np.array([1.25,1.75])*np.sin(-np.pi/8), line_color='black', 
                mode='lines')

fig.add_annotation(x=1.2, y=0, text=r'$\require{color}{\color[rgb]{0.315209,0.728565,0.037706}p}$', ax=0.5, ay=0, xref='x', 
                   yref='y', axref='x', ayref='y', arrowhead=5)
fig.add_annotation(x=1.8, y=0, text=r"$\require{color} {\color[rgb]{0.315209,0.728565,0.037706}p} + \
                   \frac{d {\color[rgb]{0.315209,0.728565,0.037706}p}}{dr} \delta r$", ax=3, ay=0, xref='x', yref='y', 
                   axref='x', ayref='y', arrowhead=5)
fig.add_annotation(x=0.6, y=0.75, ax=-0.05, ay=-0.09, xref='x', yref='y', axref='x', ayref='y',arrowhead=5)
fig.add_scatter(x=[0.25], y=[0.5], mode='text', text=r'$r_1$')
fig.add_annotation(x=0.6, y=-1.9, ax=-0.05, ay=0.09, xref='x', yref='y', axref='x', ayref='y', arrowhead=5)
fig.add_scatter(x=[0.25], y=[-0.5], mode='text', text=r'$r_2$')

fig.update_layout(yaxis=dict(scaleanchor="x", scaleratio=1, showticklabels=False), xaxis=dict(showticklabels=False),
                  showlegend=False, plot_bgcolor='rgba(0, 0, 0, 0)')
fig.show()

To find how the velocity changes with radius we need to use Newton's second law of motion: $Force = \require{color}{\color[rgb]{0.501961,0.250953,0.010028}m}a$ \ The only force acting on the flow is pressure since the flow is inviscid and therefore does not experience shear force \ The flow changes direction as it moves around the bend meaning there must be a radial acceleration and therefore a force causing the flow to move in a circular path around the origin \ This acceleration/force is centripetal and is given by $a_{c} = -\frac{{\color[rgb]{0.059472,0.501943,0.998465}v}^{2}(r)}{r} \hat{r}$ \ To find the mass we can use the fact that mass is equivalent to density times volume, when done per unit width this is $\require{color}{\color[rgb]{0.501961,0.250953,0.010028}m}={\color[rgb]{0.814433,0.253157,0.091125}\rho}(\delta l \delta r)$ \ Then we can use the force diagram above and write out Newton's second Law: $$ \require{color}{\color[rgb]{0.315209,0.728565,0.037706}p} \delta l -\left({\color[rgb]{0.315209,0.728565,0.037706}p} + \frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{dr} \delta r\right)\delta l = -{\color[rgb]{0.814433,0.253157,0.091125}\rho}(\delta l \delta r) \frac{{\color[rgb]{0.059472,0.501943,0.998465}v}^{2}(r)}{r} $$ This can then be solved for $\require{color}\frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{dr}$: $$ \require{color}\frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{dr} = {\color[rgb]{0.814433,0.253157,0.091125}\rho} \frac{{\color[rgb]{0.059472,0.501943,0.998465}v}^{2}(r)}{r} $$ From this we can see that the normal pressure gradient is what causes the flow to move in a circular path around the bend \ To solve for velocity we will need Bernoulli's equation: $\require{color}{\color[rgb]{0.315209,0.728565,0.037706}p}_{0}={\color[rgb]{0.315209,0.728565,0.037706}p}+\frac{1}{2}{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^{2}$ where $\require{color}{\color[rgb]{0.315209,0.728565,0.037706}p}_{0}$ is a constant\ In order to find an expression for the pressure gradient in terms of velocity we must take the derivative of Bernoulli's equation with respect to radius as that is the only variable that changes velocity and pressure, remembering that this is incompressible flow we get: $\require{color}0=\frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{dr}+ {\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}\frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}}{dr}$, plugging in our previous result for $\require{color}\frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{dr}$ we have: $\require{color}0={\color[rgb]{0.814433,0.253157,0.091125}\rho} \frac{{\color[rgb]{0.059472,0.501943,0.998465}v}^{2}(r)}{r} + {\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}(r) \frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}(r)}{dr}$ \ Next, we can separate the varaibles and integrate both sides: $\require{color}-\frac{dr}{r} = \frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}(r)}{{\color[rgb]{0.059472,0.501943,0.998465}v}(r)} \longrightarrow {\color[rgb]{0.059472,0.501943,0.998465}v}(r)=\frac{C}{r}$ where C is a constant to be determined by initial conditions \ From this we can see that the velocity is only dependent on the radius so velocity will be constant along at each radius, as we expected \ Using the boundary conditions that $\require{color}{\color[rgb]{0.059472,0.501943,0.998465}v}(r_1)={\color[rgb]{0.059472,0.501943,0.998465}v}_0$ we can get that $\require{color}C={\color[rgb]{0.059472,0.501943,0.998465}v}_{0}r_{1}$ therefore $\require{color}{\color[rgb]{0.059472,0.501943,0.998465}v}(r)={\color[rgb]{0.059472,0.501943,0.998465}v}_{0}\frac{r_{1}}{r}$ where $r_{1}$ is the inner radius of the bend

In [4]:
# define necessary variables to calculate the speed at each radii
v_0 = 100 # initial speed at r_1
r_1 = 1 # inner radius
r_2 = 2 # outer radius
In [5]:
# calculate the speed between r_1 and r_2
r = np.linspace(r_1, r_2, 9)
v = v_0*r_1*(1/r)

To solve for pressure we must use Bernoulli's equation again\ If we assume that at $r_1$ the pressure is atmospheric pressure $\require{color}{\color[rgb]{0.315209,0.728565,0.037706}p}_a = 101.325 \ kPa$ then we can find the stagnation pressure which will be constant at any point in the flow using Bernoulli's equation $\require{color}{\color[rgb]{0.315209,0.728565,0.037706}p}_{0}={\color[rgb]{0.315209,0.728565,0.037706}p}+\frac{1}{2}{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^{2}$ and we can find the dynamic pressure term from how we solved for velocity giving: $\require{color}{\color[rgb]{0.315209,0.728565,0.037706}p}_0 = {\color[rgb]{0.315209,0.728565,0.037706}p}_{\color[rgb]{0.989013,0.435749,0.811750}a} + \frac{1}{2}{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}_0^2$ $$ \require{color}{\color[rgb]{0.315209,0.728565,0.037706}p} = {\color[rgb]{0.315209,0.728565,0.037706}p}_0 - \frac{1}{2}{\color[rgb]{0.814433,0.253157,0.091125}\rho} \frac{C^2}{r^2} = {\color[rgb]{0.315209,0.728565,0.037706}p}_0 - \frac{1}{2}{\color[rgb]{0.814433,0.253157,0.091125}\rho} \ \frac{{\color[rgb]{0.059472,0.501943,0.998465}v}_0^2r_1^2}{r^2} \longrightarrow {\color[rgb]{0.315209,0.728565,0.037706}p} = {\color[rgb]{0.315209,0.728565,0.037706}p}_{\color[rgb]{0.989013,0.435749,0.811750}a} + \frac{1}{2} {\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}_0^2\left(1-\frac{r_1^2}{r^2}\right) $$ This equation will find the pressure at any radius

In [6]:
# define variables needed to calculate pressure
rho = 1.293 # denisty of air in kg/m^3
p_a = 101325 # pressure of air (atmospheric pressure) in Pa
p_0 = p_a + (1/2)*rho*v_0**2
In [7]:
# calculate pressure between r_1 and r_2
p = p_0 - (1/2)*rho*(v_0**2)*(r_1**2)/(r**2)
In [8]:
# define variables needed to make a contour plot
theta = np.linspace(0, np.pi/2, 50)
x=r_2*np.cos(theta)
y=np.linspace(-2,2,50)
In [9]:
# define pressure and velocity for contour plot
v_z = []
p_z = []
for y_val in y:
    v_y_list = []
    p_y_list = []
    for x_val in x:
        if (x_val**2 + y_val**2) < r_1:
            p_xy = p_a
            v_xy = v_0
        else:
            v_xy = v_0/np.sqrt(x_val**2 + y_val**2)
            p_xy = p_0 - (1/2)*rho*(v_0**2)*(r_1**2)/(x_val**2 + y_val**2)
        v_y_list.append(v_xy)
        p_y_list.append(p_xy)
    v_z.append(v_y_list)
    p_z.append(p_y_list)
In [10]:
# make contour plots for velocity and pressure  
fig = go.Figure()
fig.add_contour(z=v_z, x=x, y=y, contours_coloring='heatmap', line_width=0, colorbar=dict(title='Velocity (m/s)', 
                titleside='right'), colorscale='PuBuGn', reversescale=True)
fig.add_contour(z=p_z, x=x, y=y, contours_coloring='heatmap', line_width=0, visible=False, colorbar=dict(title=
                'Pressure (Pa)', titleside='right'), colorscale='dense', reversescale=True)

fig.add_shape(type='circle', xref='x', yref='y', x0=-1, y0=-1, x1=1, y1=1, fillcolor='white', line_color='white')
fig.add_scatter(x=2*np.cos(t2), y=2*np.sin(t2), line_width=0)
fig.add_scatter(x=3*np.cos(t2), y=3*np.sin(t2), line_color='white', fill='tonexty', fillcolor='white')

fig.update_layout(yaxis=dict(scaleanchor="x", scaleratio=1, showticklabels=False, range=[-2,2]), 
                  xaxis=dict(showticklabels=False), plot_bgcolor='rgba(0, 0, 0, 0)', autosize=False, width=800, height=800, 
                  showlegend=False, updatemenus=[dict(active=0, buttons=[dict(label='Velocity', method='update', 
                  args=[{'visible':[True, False, True, True, True]}]), dict(label='Pressure', method='update',
                  args=[{'visible':[False, True, True, True, True]}])], y=1.2, x=0.41)])
fig.show()

As you can see in the contour plots above the pressure and velocity are inversely proportional to one another which is what we would expect \ Additionally, the velocity decreases with increasing radius which we would also expect because there is less centripetal acceleration acting on the flow the farther you get from the center of the bend

In [11]:
# define lines where we want to visualize the flow and time and distance information for those lines
total_distance = np.pi*r
total_time = total_distance / v
dist = total_time[7] * v

r_lines = np.array([r[1], r[4], r[7]])
r_dist = np.array([dist[1], dist[4], dist[7]])
rad = -np.pi/2 + (r_dist / r_lines)
In [12]:
# create animated plot to show ideallized flow movement 
fig = go.Figure()

frames = [go.Frame(data=[]) for k in range(50)]
for i, r_line in enumerate(r_lines):
    fig.add_scatter(x=r_line*np.cos(t2), y=r_line*np.sin(t2), line_dash='dash', line_color='lightsteelblue')
    x_pos = r_line*np.cos(np.linspace(-np.pi/2, rad[i], 50))
    y_pos = r_line*np.sin(np.linspace(-np.pi/2, rad[i], 50))
    try:
        end = np.where(x_pos < 0)[0][0]
        x_pos[end:len(x_pos)] = 'NaN'
    except IndexError:
        None
    for j, f in enumerate(frames):
        f.data += (go.Scatter(x=[x_pos[j]], y=[y_pos[j]], mode='markers', marker=dict(color='teal', size=10)),)

fig.update(frames=frames)
fig.update_layout(updatemenus = [dict(type='buttons', buttons=[dict(label='play', method='animate', 
                  args=[None, {'frame':{'duration':100}}])])], xaxis=dict(range=[0,2], showticklabels=False), 
                  yaxis=dict(scaleanchor="x", scaleratio=1, showticklabels=False), showlegend=False, 
                  plot_bgcolor='rgba(0, 0, 0, 0)', autosize=False, width=900, height=700)

fig.add_scatter(x=2*np.cos(t2), y=2*np.sin(t2), fill='tozeroy', fillcolor='rgba(133, 163, 201, 0.5)', line_color='black')
fig.add_scatter(x=np.cos(t2), y=np.sin(t2), fill='tozeroy', fillcolor='white', line_color='black')

fig.show()
In [ ]: